{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nTime advection\n==============\n\nExample which use CMEMS surface current with a Runge-Kutta 4 algorithm to advect particles.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 2\nimport re\nfrom datetime import datetime, timedelta\n\nfrom matplotlib import pyplot as plt\nfrom matplotlib.animation import FuncAnimation\nfrom numpy import arange, isnan, meshgrid, ones\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.dataset.grid import GridCollection\nfrom py_eddy_tracker.gui import GUI_AXES\n\nstart_logger().setLevel(\"ERROR\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class VideoAnimation(FuncAnimation):\n    def _repr_html_(self, *args, **kwargs):\n        \"\"\"To get video in html and have a player\"\"\"\n        content = self.to_html5_video()\n        return re.sub(\n            r'width=\"[0-9]*\"\\sheight=\"[0-9]*\"', 'width=\"100%\" height=\"100%\"', content\n        )\n\n    def save(self, *args, **kwargs):\n        if args[0].endswith(\"gif\"):\n            # In this case gif is used to create thumbnail which is not used but consume same time than video\n            # So we create an empty file, to save time\n            with open(args[0], \"w\") as _:\n                pass\n            return\n        return super().save(*args, **kwargs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Data\n----\nLoad Input time grid ADT\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "c = GridCollection.from_netcdf_cube(\n    get_demo_path(\"dt_med_allsat_phy_l4_2005T2.nc\"),\n    \"longitude\",\n    \"latitude\",\n    \"time\",\n    # To create U/V variable\n    heigth=\"adt\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Anim\n----\nParticles setup\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "step_p = 1 / 8\nx, y = meshgrid(arange(13, 36, step_p), arange(28, 40, step_p))\nx, y = x.reshape(-1), y.reshape(-1)\n# Remove all original position that we can't advect at first place\nt0 = 20181\nm = ~isnan(c[t0].interp(\"u\", x, y))\nx0, y0 = x[m], y[m]\nx, y = x0.copy(), y0.copy()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Function\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def anim_ax(**kw):\n    fig = plt.figure(figsize=(10, 5), dpi=55)\n    axes = fig.add_axes([0, 0, 1, 1], projection=GUI_AXES)\n    axes.set_xlim(19, 30), axes.set_ylim(31, 36.5), axes.grid()\n    line = axes.plot([], [], \"k\", **kw)[0]\n    return fig, axes.text(21, 32.1, \"\"), line\n\n\ndef update(_):\n    tt, xt, yt = f.__next__()\n    mappable.set_data(xt, yt)\n    d = timedelta(tt / 86400.0) + datetime(1950, 1, 1)\n    txt.set_text(f\"{d:%Y/%m/%d-%H}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "f = c.filament(x, y, \"u\", \"v\", t_init=t0, nb_step=2, time_step=21600, filament_size=3)\nfig, txt, mappable = anim_ax(lw=0.5)\nani = VideoAnimation(fig, update, frames=arange(160), interval=100)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Particules stat\n---------------\nTime_step settings\n^^^^^^^^^^^^^^^^^^\nDummy experiment to test advection precision, we run particles 50 days forward and backward with different time step\nand we measure distance between new positions and original positions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure()\nax = fig.add_subplot(111)\nkw = dict(\n    bins=arange(0, 50, 0.002),\n    cumulative=True,\n    weights=ones(x0.shape) / x0.shape[0] * 100.0,\n    histtype=\"step\",\n)\nkw_p = dict(u_name=\"u\", v_name=\"v\", nb_step=1)\nfor time_step in (10800, 21600, 43200, 86400):\n    x, y = x0.copy(), y0.copy()\n    nb = int(30 * 86400 / time_step)\n    # Go forward\n    p = c.advect(x, y, time_step=time_step, t_init=20181.5, **kw_p)\n    for i in range(nb):\n        t_, _, _ = p.__next__()\n    # Go backward\n    p = c.advect(x, y, time_step=time_step, backward=True, t_init=t_ / 86400.0, **kw_p)\n    for i in range(nb):\n        t_, _, _ = p.__next__()\n    d = ((x - x0) ** 2 + (y - y0) ** 2) ** 0.5\n    ax.hist(d, **kw, label=f\"{86400. / time_step:.0f} time step by day\")\nax.set_xlim(0, 0.25), ax.set_ylim(0, 100), ax.legend(loc=\"lower right\"), ax.grid()\nax.set_title(\"Distance after 50 days forward and 50 days backward\")\nax.set_xlabel(\"Distance between original position and final position (in degrees)\")\n_ = ax.set_ylabel(\"Percent of particles with distance lesser than\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Time duration\n^^^^^^^^^^^^^\nWe keep same time_step but change time duration\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure()\nax = fig.add_subplot(111)\ntime_step = 10800\nfor duration in (10, 40, 80):\n    x, y = x0.copy(), y0.copy()\n    nb = int(duration * 86400 / time_step)\n    # Go forward\n    p = c.advect(x, y, time_step=time_step, t_init=20181.5, **kw_p)\n    for i in range(nb):\n        t_, _, _ = p.__next__()\n    # Go backward\n    p = c.advect(x, y, time_step=time_step, backward=True, t_init=t_ / 86400.0, **kw_p)\n    for i in range(nb):\n        t_, _, _ = p.__next__()\n    d = ((x - x0) ** 2 + (y - y0) ** 2) ** 0.5\n    ax.hist(d, **kw, label=f\"Time duration {duration} days\")\nax.set_xlim(0, 0.25), ax.set_ylim(0, 100), ax.legend(loc=\"lower right\"), ax.grid()\nax.set_title(\n    \"Distance after N days forward and N days backward\\nwith a time step of 1/8 days\"\n)\nax.set_xlabel(\"Distance between original position and final position (in degrees)\")\n_ = ax.set_ylabel(\"Percent of particles with distance lesser than \")"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}